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Abstract 

We generated configurations with the parametrized fixed-point Dirac 
operator Dpp on a (1.6 fm) 4 box at a lattice spacing a — 0.13 fm. We 
compare the distributions of the three lowest k = 1, 2, 3 eigenvalues in the 
v — 0, 1, 2 topological sectors with that of the Random Matrix Theory 
predictions. The ratios of expectation values of the lowest eigenvalues 
and the cumulative eigenvalue distributions are studied for all combina- 
tions of k and v. After including the finite size correction from one-loop 
chiral perturbation theory we obtained for the chiral condensate in the 
MS scheme [E(2GeV)] 1/3 = 0.239(11) GeV, where the error is statistical 
only. 



Spontaneous chiral symmetry breaking and the related existence of light 
Goldstone bosons is a basic feature of QCD. Chiral Perturbation Theory (ChPT) 
provides a systematic description of this physics in terms of a set of low en- 
ergy constants which encode the related non-perturbative features of QCD. The 
method of low energy effective Lagrangians simplifies the calculations signifi- 
cantly [H |2l El H] and over the years ChPT became a refined powerful technique. 
Gasser and Leutwyler recognized very early, decades before the numerical cal- 
culations could attack such problems, that these constants can also be fixed 
using physical quantities which, presumably, will never be measured in real 
experiments. They can be studied, however, in lattice QCD. 

The e- regime [5j [6l [3 [8] describes physics close to the chiral limit in a box 
whose size is larger than the QCD scale. On the other hand, the size of the 
box relative to the Goldstone boson correlation length must be small. Under 
these conditions the Goldstone bosons, as opposed to other excitations, feel 
the effect of boundaries strongly. ChPT provides a powerful and systematic 
way to calculate the finite size corrections. A nice additional feature of the e- 
regime is that Random Matrix Theory (RMT) |9J makes precise predictions for 
microscopic observables. RMT relates in particular the distribution of low-lying 
eigenvalues of the Dirac operator in different topological sectors to the chiral 
condensate. 

The pioneering numerical works [TQl [TTJ [12l [13l [14j [15l US] in the e-regime 
in quenched QCD suggested that this regime can be an excellent tool to study 
low-energy physics in QCD. The special problems of the e-regime called for 
new numerical procedures [HI [15l HH [17] which were first tested also in this 
approximation. Combining these numerical developments with the renormal- 
ization properties of the spectral density the work [18j present a state of the 
art analysis for the quark condensate and the first Leutwyler-Smilga sum rule 
in quenched QCD. The quenched approximation, however, tends to be singular 
in the chiral limit and it is expected that quenching is even more problematic 
here than in other cases [19] . 

Unfortunately, full QCD simulations are expensive. In addition, in the 
e— regime the Dirac operator should have excellent chiral properties. The stan- 
dard choice is the overlap Dirac operator [20] with a hybrid Monte Carlo algo- 
rithm. 

The results in [21] reflect already the basic physics features of the e-regime. 
The distribution of the low-lying eigenmodes could be fitted quite well to the 
RMT predictions and a reasonable value for the chiral condensate was obtained. 
These results are promising given the fact that the box was small (1.3 fm), 
the quark mass was rather large (m q > 40MeV) and the lattice was coarse 
(a = 0.16 fm). The results in a larger box of 1.5 fm and at smaller quark mass 
m q 20 MeV obtained in [22] were consistent with those of [21] . 
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The first serious simulation results in the e-regime have been presented by the 
JLQCD group recently [23] ■ In that work Nf = 2 QCD was simulated with over- 
lap fermions (created with the Wilson kernel) on a lattice of size 1.8 3 x 3.6 fm 4 
with a resolution a = 0.11 fm obtained from rQ. Different bare quark masses 
(am q ) were considered in the range 0.110, . . . , 0.020, 0.002. The smallest quark 
mass (corresponding to w 3MeV) was certainly small enough to reach the e- 
regime. Hybrid Monte Carlo with overlap fermions has problems when the topo- 
logical charge changes. To avoid this an action was used which prevented topol- 
ogy change and the whole run stayed in the Q = sector. The authors observed 
an overall good agreement with RMT. The fermion condensate, which is the 
only parameter to be fitted, was found to be E(2GeV) 1 / 3 = 0.251(7)(ll)GeV 
in the MS scheme. 

Our work is based on the parametrized fixed-point (FP) action which has 
been tested in detail in quenched QCD [24]. The exact FP Dirac operator 
satisfies the Ginsparg- Wilson relation 

4w + d gw = Dl w 2RD GW , (1) 

where R is a local operator and is trivial in Dirac space. The parametrized FP 
action has many gauge paths and involves a special smearing with projection 
to the gauge group SU(3). As a consequence, hybrid Monte Carlo algorithms 
can not be used. In the work [55] we advised a partially global update with 
three nested accept/reject steps which can reach small quark masses even on 
coarse lattices. Several steps of this algorithm were developed in [26] using 
some earlier suggestions (27] [28] |29j [30] . All the three steps are preconditioned. 
Pieces of the quark determinant are switched on gradually in the order of their 
computational expenses. The largest and most fluctuating part of the deter- 
minant is coming from the ultraviolet modes. We reduce these fluctuations by 
calculating the trace of D% P ,n = 1,...,4 [3U [32 [26] . The ~ 100 lowest-lying 
modes are calculated and subtracted. The determinant of the reduced and sub- 
tracted Dpp is calculated stochastically [33] [331 [351 [35] . For the determinant 
breakup we generalized the mass shifting method of [3T] • For the strange quark 
we performed a root operation. Since the strange quark mass is not very small, 
we used a polynomial expansion to approximate this root operator. Recently 
we replaced the polynomial expansion by a rational approximation 37J which 
brought a rs 30% performance gain in this part. 

Due to the partial global updating procedure the algorithm, beyond a cer- 
tain volume, scales with V 2 which constrains the size of lattices which can be 
considered. Due to the large number of subtracted low-lying modes, small quark 
mass is not a barrier. 

We generated ~ 4000 configurations in the Markov chain on an 12 4 lattice us- 
ing the partially global algorithm discussed above. The distance between every 
second configurations in this chain is similar to that of two gauge configurations 
separated by a typical Metropolis gauge update sweep. We considered every 
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tenth configurations from the Markov chain and performed the measurements 
on the remaining ~ 400 configurations. When calculating statistical errors we 
formed 20 bins and used the jackknife method. 



We fixed the lattice spacing from the Sommer parameter 
and found a — 0.129(5) fm. Our box has the size (1.6 fm) 4 . 



r = 0.49 fm 



The Dirac operator -Dfp has no exact chiral symmetry due to parametriza- 
tion errors. As a consequence, the quark masses have an additive mass renormal- 
ization. The degenerate u, d and the s quark masses in the code are M uc j = 0.025 
and M a — 0.103. The additive renormalization was measured following the steps 
described in [39] and has the value M = 0.0147(3), as shown in Fig. [TJ Sub- 
tracting the additive renormalization we get the bare masses m u( j = 0.0103(3) 
and m s = 0.0883(3) which corresponds to 16MeV and 137 MeV, respectively. 



0.08 




Figure 1: The AWI mass for the three combinations of quark masses. The 
linear extrapolation to Mawi — gives an additive mass renormalization Mq = 
0.0147(3). 



As we mentioned above, a simulation with lighter quark masses would not 
be more expensive. In our case, however, a smaller to u( j were not useful cither. 
RMT predicts the probability distribution p v k(£uk) for the /c-th low- lying eigen- 
value (k — 1, 2, . . . ) of the Dirac operator in the topological sector v. Denoting 
the corresponding eigenvalues of the (continuum) Dirac operator by ia u k , the 
variable is related to the bare chiral condensate X as = a„fcXV ". Here V 
is the volume of the box. Fig.[2]shows the prediction of RMT for the cumulative 
distributions f^" 1 * d£jVfc(£) for v — and k = 1, 2, 3. The distributions depend 
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on /Ltj = nu'EV where are the quark masses. Decreasing the to u( j mass by 
more than a factor 10 (at fixed m s ) the cumulative distributions practically 
remain unchanged. On the other hand sending the strange quark mass to in- 
finity (at fixed m u d) we land on a Nf = 2 flavor theory with a visibly different 
distribution. The strange quark has a (modest) effect on our observables. 




Figure 2: Random Matrix Theory prediction for the cumulative distribution 
of £ = = ctykYjV , where ia v k is the k-th eigenvalue of the continuum Dirac 
operator in a gauge background with topological charge v, and m = rriiEV. 
Here the v = results are shown, but the picture is similar for u = 1,2. 



While the eigenvalues of the continuum Dirac operator lie on the imagi- 
nary axis the spectrum of lattice Dirac operators is more complicated. As it 
is well known, the spectrum of the Dirac operator satisfying the Ginsparg- 
Wilson relation with 2R = 1 lies on the circle |A — 1| = 1. In our case 
where 2R is different from 1, it is convenient to introduce the rescaled oper- 
ator -Dgw = V2i?Z?GW\/2i?, for which the 2R factor is eliminated in the GW 
relation, and the spectrum lies on the circle. (In fact, our operator 2R for the 
low- lying modes is effectively a constant close to 1 within a few percent.) To 
relate the eigenvalues on the GW circle to those appearing in the RMT (or in 
general, in continuum expressions) it is natural to use the stereographic projec- 
tion 

ia= i^xT2- < 2 > 

As Fig. [3] shows, the low- lying eigenvalues of Dpp( m ) (which we obtain during 
our simulation) are close but not exactly on the GW circle, hence we project 



5 



them first horizontally onto the GW circle, and make the stereographic pro- 
jection in the next step. This is justified by the observation that making a 
systematic GW improvement of -Dfp towards Dqw the imaginary part of the 
eigenvalues stays practically constant - they move horizontally to the GW cir- 
cle. The values of a obtained by this procedure are then compared with the 
RMT predictions. 
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Figure 3: The histogram of the deviation from the circle, 6r = 1 — |A — 1|, for 
the eigenvalues of -Dfp on 8 3 x 24 and 12 4 lattices for the complex eigenvalues 
with |ImA| < 0.1, using the present mass parameters. 



Having the probability distributions from RMT we can calculate the ratios 
(£uk)/(£u'k'), where the factor Y,V cancels. These predictions are compared 
with the measured ratios in Fig.[3J where all the ratios in the topological sectors 
v = 1, 2, 3 for the first three lowest eigenvalues are shown. 

At this point we have to discuss the way we identify the topological sectors. 
According to the index theorem [40j one can identify the topological charge from 
the zero modes of a Ginsparg- Wilson Dirac operator. However, our Dfp has no 
exact chiral symmetry and has modes on the real axis. Occasionally, the real 
eigenvalue A might be even far from zero in which case the topological interpre- 
tation is uncertain. We note here that using an exact Ginsparg- Wilson operator 
(with kernel -Dfp) for measuring the topological charge is not a good solution 
to this problem. The overlap just projects most of the real eigenvalues to the 
point A — 0. Following the intuitive picture that topology is related to extended 
objects we investigated the correlation between the inverse participation ratio 
(IPR) of the eigenvector having a real eigenvalue A and the value of A. Here 
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Figure 4: Ratios of expectation values (a„t) obtained from simulations com- 
pared to the results of RMT, where vk is indexing the fc-th lowest eigenvalue in 
the topological sector v. The different symbols refer to the denominator while 
vk of the numerator is indicated at the data points. For example, the highest 
ratio with value ~ 3 in the figure refers to (^23) /(^oi) 



IPR is given by Y^l=i l^i ( x )\ 4 f° r a normalized eigenvector ip( x \ The IPR 
is inversely proportional to the size of the effective support of the eigenfunction. 
As Fig. [5] shows, there is a strong correlation indeed. As A is moving away 
from zero, the IPR increases indicating that the wave function becomes more 
and more localized. Fig. [6] demonstrates that the real eigenvalues A(Z?fp) are 
strongly concentrated in a small region close to zero. The histogram of A(Dfp) 
has a large, narrow peak followed by a long tail. 

Fig. [7| shows the histogram of IPR for the real cigcnmodes of Z?fp- There 
is a strong, narrow peak which corresponds to extended real modes and a long 
tail corresponding to wave functions with smaller support. 

For the questions we study in this paper it is not necessary to use a procedure 
which leaves the relative weights of the different topological sectors intact. (Due 
to the long autocorrelation time of the topological charge it would be anyhow 
difficult to control these weights.) We introduce therefore a cut above which 
the real eigenvalues are discarded. This cut must be sufficiently small to allow 
for a good identification of the topological charge, but large enough to allow for 
a reasonable number of configurations in each sectors we want to investigate. 
The cut should be smaller than the spectral gap (for the spectrum see Fig. [8]), 



7 



0.004 



o I ^V--? - i , i , i i i , 

0.1 0.2 0.3 0.4 0.5 

Figure 5: The inverse participation ratio y^, 1 _ 1 O^)] 4 f° r reaZ eigenvalues 
A(£>fp) vs. the eigenvalue. Eigenvectors with larger real part have smaller 
extension in agreement with the expectations. 



otherwise the effect of topology will be strongly reduced. We introduced a cut of 
size 0.03 which is approximately equal to the minimal gap in the v = 1 sector. 
We accepted a real eigenvalue A as a sign of topology if A < 0.03. This cut 
splits our sample of 414 configurations into 328, 51 and 35 configurations for 
the v = 0, 1, 2 sectors, respectively. This cut corresponds also to the end of the 
peak of the IPR distribution shown in Fig. [7J 

Fig. [U compares the cumulative distributions of the first three eigenvalues 
in the v = topological sector obtained using the overlap improved fixed-point 
Dirac operator -Dqw vs - ^fp m the measurement. The distributions are very 
similar showing that the two operators are close to each other. The k = 1 
distribution from -Dqw has a tail at small eigenvalues (i.e. has a smaller gap than 
-Dfp) demonstrating a systematic error one obtains when using different Dirac 
operators for the generation of the configurations and for the measurements: 
the small eigenvalues are less effectively suppressed. 

In Figs. [TUl HU HU the cumulative distributions of the D-pp operator are 
compared with RMT for the v = 0, 1, 2 topological sectors. The only matching 
parameter is the bare condensate E which enters the RMT predictions through 
/ij = rriiY^V and in aY^V. We obtained the result for this bare quantity E 1 / 3 = 
0.291(3)(9) GeV. Both errors are statistical. The first error comes from the 
statistical errors of the measured distributions, the second one is due to the 
error in the lattice scale a. 
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Figure 6: Histogram of A(Z?fp) for the real eigenvalues. In determining the 
topological charge v we considered only the real eigenvalues with A < 0.03. 



The result above receives a small correction due to the fact that, for simpli- 
fying the presentation, we suppressed a technical complication. As mentioned 
before the exact fixed-point operator satisfies a Ginsparg- Wilson relation with 
a local operator 2R. For this reason the quark mass enters in a simple additive 
way in the form D + m(l - D/2), where D = V2RD FP V2R. Effectively the 
operator 2R behaves in the infrared like a constant close to 1. Its expectation 
value for the lowest ~ 100 eigenvectors is 1.05 within 1%. Using the spectrum 
of D (which is in fact identical to that of 2RDpp) the matching with RMT gives 
the slightly changed result E 1 / 3 = 0.286(3) (9) GeV. 

In the language of a ferromagnetic O(N) model one should interpret S as 
the absolute value of the magnetization in the finite volume V. This differs from 
the value £oo defined in the infinite volume by a finite-size correction which can 
be calculated in ChPT. In the presence of the magnetic field h the orientation 
of magnetization is controlled by the Boltzmann factor exp(hM cos 8) where M 
is the total magnetization. This gives 

where Yn(z) is related to the modified Bessel functions. Comparing this ex- 
pression with the result of ChPT [7] we get for TV = 4 (corresponding to two 
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Figure 7: Histogram of inverse participation ratio for the real eigenmodes of 
Dpp. 



flavors) E = pEoc, where 



Here the shape coefficient /3i takes the value 0.14046 for a symmetric box [7J. 
The one-loop finite-size correction to the order parameter (magnetization in the 
example above) has been calculated for the SU(iV/) x SU(iV/) symmetry group 
in [51 [B], for the O(N) group in [IT] and up to the two-loop level in [7J. 

Neglecting the strange quark contribution in the finite size effects we use the 
two-flavor result to correct the measur ed bare E 1 / 3 = 0.286(3) (9) GeV by the 
one-loop finite size correction to obtain E^ 3 = E 1 / 3 /1.119 = 0.255(3)(9) GeV. 

In the previous version of this paper we quoted only this bare value. Since 
then we have calculated the renormalization factor connecting the bare lattice 
result to the MS scheme. The main points of this calculation are summarized 
in the Appendix. With the conversion factor 0.82(3) we obtain 




(4) 



[Sj^(2 GeV)] 1/3 = 0.239(11) GeV. 



(5) 
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Figure 8: The spectrum of Dpp for 50 configurations on 8 3 x 24 and 12 lattices, 
at the same bare mass parameters and gauge coupling. 




Figure 9: The cumulative distributions in the v = topological sector for the 
three lowest complex eigenvalues as measured with Z?fp and with the overlap 
operator Dgw using the fixed-point operator as kernel. 
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Figure 10: The cumulative distribution of = a^tSl/ for v = 0, fc = 1, 2, 3. 
Here a is obtained from A(-Dfp) by stereographic projection described in the 
text, and is rescaled using S 1//3 = 0.291(3) GeV. The continuous line is the 
RMT prediction at /Mi = f^2 = m u dT,V = 1.43 and /13 = m s 'EV = 12.3. 




Figure 11: The cumulative distribution of £„k = a„fcEV for v = 1, k = 1, 2, 3. 
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APPENDIX 



The bare lattice scalar density is not universal, it depends on the lattice action 
and it needs renormalization. The conventional way is to express the result in 
MS scheme at some given scale, e.g. fio = 2 GeV. 

To relate the bare lattice result to MS scheme one generally uses the Rome- 
Southhampton RI/MOM method [42]. In this method one introduces an inter- 
mediate scheme RI (or RI0) where some Green's functions of the actual operator 
(calculated in a fixed gauge) are equated to their corresponding Born terms at 
a given scale p 2 = /i 2 . 

In the first step one relates the bare lattice operator non-perturbatively to 
its renormalized counterpart in the RI scheme 

OgM=Z^MO& e {a). (A-l) 

The matching scale ^ is restricted by two conditions: It should be sufficiently 
large to avoid non-perturbative effects and small enough to avoid large cut-off 
effects. For the scalar density, the non-perturbative matching was done in [33], 
using a technique which removes a large part of 0(a) cut-off effects. 

In the second step one relates the renormalized operators in the RI and MS 
schemes in perturbation theory 

0fW = ^ RI «M. (A-2) 

Obviously, the value of /i must lie in the perturbative regime. The factor 
Z m ' m (fx) has been calculated up to NNLO [44 and NNNLO [45]. Combin- 
ing (|A-lj) and (|A-2|) . we obtain the matching factor connecting the bare lattice 
and the renormalized MS results 

Z^ lat ( M , a) = Z^^MZg^fji, a) . (A-3) 



As the expansion parameter in the perturbative series is the running coupling 
Qt(fJk) of QCD, we need to determine it for the scale \x where the matching is 

1 These schemes differ in their definition of the quark field renormalization factor Z q . 
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done. The running of a is described by the differential equation 



where the 4-loop (3- function in MS scheme is given in [48]. As initial condition 
we use a(2 GeV) = 0.2904 from [49] @ .The procedure described above is carried 
out at several matching scales fi > 2 GeV. 

In the third step Renormalization Group technique is used to obtain the 
desired conversion factor Z s MS ' lat (/io, a) for the scalar density at /zq = 2 GeV, 



£ s MS ' lat ( M o,a) = Z s MS >Va)exp ( - f " da^ ) . (A-5) 

V -Wo) /3(a) J 

Here, due to the relation Z s = , the anomalous dimension j s is related to 
the mass anomalous dimension by 7 S = — j m . The latter has been calculated to 
four loops in [H] g7] . 

In Fig. [T3] we plot our result Z^ S ' lat (/xo , o) at /io = 2 GeV vs. the matching 
scale fj, for the intermediate schemes RI and RI . If the non-perturbative, cut- 
off and higher order perturbative effects were negligible, the two curves would 
coincide and would be independent of the choice of fi. Taking the average values 
for 2 GeV < [i < 3 GeV and taking into account the slight difference between 
the curves, we get 

Zf^'^ifio = 2 GeV, a) = 0.82(3) (A-6) 



2 This result is obtained by using the PDG value a(M z ) = 0.1176 ± 0.002 [50] and running 
it from 5-fiavors down to 3-flavors, across the and m c thresholds 
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Figure 13: The renormalization factor Z s MS ' lat (/io, a) for the scalar density, 
connecting the MS and lattice schemes, at /j,q = 2 GeV, obtained using RI and 
RI' schemes at the intermediate step vs. the matching scale \i. 



16 



References 



[1] S. Weinberg, Physica A96 (1979) 327. 

[2] J. Gasser and H. Leutwyler, Phys. Lett. B 125 (1983) 321, 325. 

[3] J. Gasser and H. Leutwyler, Ann. Phys. (N.Y.). 158 (1984) 142. 

[4] J. Gasser and H. Leutwyler, Nucl. Phys. B 250 (1985) 465. 

[5] J. Gasser and H. Leutwyler, Phys. Lett. B 184 (1987) 83; 188 (1987) 477. 

[6] J. Gasser and H. Leutwyler, Nucl. Phys. B 307 (1988) 763. 

[7] P. Hasenfratz and H. Leutwyler, Nucl. Phys. B 343 (1990) 241. 

[8] F. C. Hansen, Nucl. Phys. B 345 (1990) 685. 

[9] E. V. Shuryak and J. J. M. Vcrbaarschot, Nucl. Phys. A 560 (1993) 306 
[arXiv:hep-th/9212088], 

J. J. M. Verbaarschot and T. Wettig, Ann. Rev. Nucl. Part. Sci 50 (2000) 
343 [arXiv:hep-ph/0003017], 

P. H. Damgaard and S. M. Nishigaki, Phys. Rev. 63 (2001) 045012 
[arXiv:hep-th/0006111] 

[10] R. G. Edwards, U. M. Heller, J. E. Kiskis and R. Narayanan, 
Phys. Rev. Lett. 82 (1999) 4188. 

[11] W. Bictcnholz, K. Janscn and S. Shchcrcdin, JHEP 0307 (2003) 033. 

[12] W. Bietenholz, Th. Chiarappa, K. Jansen Kei-ichi Nagai and S. Shcheredin, 
JHEP 0402 (2004) 023. 

[13] L. Giusti, M. Liischcr, P. Wcisz and H. Wittig, JHEP 0311 (2003) 23 
[arXiv:hep-lat/0309189]. 

[14] L. Giusti, P. Hernandez, M. Laine, P. Weisz and H. Wittig, JHEP 0404 
(2004) 013 [arXiv:hep-lat/0402002]. 

[15] R. G. Edwards, Nucl. Phys. Proc. Suppl. 106 (2002) 38 [arXiv:hcp- 
lat/0111009]. 

[16] T. DeGrand and S. Schacfcr, Comput. Phys. Commun. 159 (2004) 185 
[arXiv:hep-lat/0401011]. 

[17] L. Giusti, C. Hoelbling, M. Liischer and H. Wittig, Comput. Phys. Com- 
mun. 153 (2003) 31 [arXiv:hep-lat/0212012]. 

[18] L. Giusti and S. Necco, JHEP 0704 (2007) 090 [arXiv:hcp-lat/0702013]. 

[19] P. H. Damgaard, Nucl. Phys. Proc. Suppl. 128 (2004) 47 [arXiv:hep- 
lat/0310037]. 



17 



[20] H. Neuberger, Phys. Lett. B 417 (1998) 141 [arXiv:hep-lat/9707022]; 427 
(1998) 353 [arXiv:hep-lat/9801031]. 

[21] T. DcGrand and S. Schaefer, Phys. Rev. D 71 (2005) 034507 [arXiv:hep- 
lat/0412005]; 72 (2005) 054503 [arXiv:hep-lat/0506021]; PoS LAT2005 
(2006) 140 [arXiv:hcp-lat/0508025]. 

[22] Th. DeGrand, Zh. Liu, S. Schacfcr, Phys. Rev. D 74 (2006) 094504, 
Erratum-ibid. D 74 (2006) 099904 [arXiv:hep-lat/0608019]; 

[23] H. Fukaya et al. [JLQCD Collaboration], Phys. Rev. Lett. 98, 172001 (2007) 
[arXiv:hcp-lat /0702003] . 

[24] P. Hasenfratz, S. Hauswirth, K. Holland, T. Jorg, F. Nicdermayer and 
U. Wenger, Int. J. Mod. Phys. C 12 (2001) 691 [arXiv:hep-lat/0003013], 
P. Hasenfratz, S. Hauswirth, T. Jorg, F. Niedermayer and K. Holland, 
Nucl. Phys.B 643 (2002) 280 [arXiv:hep-lat/0205010], 
C. Gattringcr et al. Bern-Graz-Regensburg Collaboration, Nucl. Phys.B 
677 (2004) 3 [arXiv:hep-lat/0307013], 

P. Hasenfratz, K. J. Juge and F. Niedermayer, Bern-Graz-Regensburg Col- 
laboration, JHEP 0412:030,2004 [arXiv:hep-lat/0411034]. 

[25] A. Hasenfratz, P. Hasenfratz and F. Niedermayer, Phys.Rev. D72 (2005) 
114508 [arXiv:hep-lat/0506024]. 

[26] A. Hasenfratz and F. Kncchtli, Comput. Phys. Commun. 148 (2002) 81 
[arXiv:hcp-lat /0203010] . 

[27] M. Grady, Phys. Rev. D 32 (1985) 1496 [arXiv:hep-lat/0203010]. 

[28] M. Creutz, Algorithms for Simulating Fermions (1992), in Quantum Fields 
on the Computer, World Scientific Publishing. 

[29] A. Alexandru and A. Hasenfratz, Phys. Rev. D 65 (2002) 114506 
[arXiv:hcp-lat/0203026]; 66 (2002) 094502 [arXiv:hcp-lat/0207014]. 

[30] F. Kncchtli and U. Wolff, (Alpha Coll.) Nucl. Phys. B 663 (2003) 3 
[arXiv:hep-lat /0303001] . 

[31] M. Hasenbusch, Phys. Lett. B 519 (2001) 177 [arXiv:hcp-lat/0107019]. 

[32] P. dc Forcrand, Nucl. Phys.proc. Suppl. 119 (1999) [arXiv:hep- 
lat/9809145]. 

[33] A. D. Kennedy and J. Kuti, Phys. Rev. Lett. 54 (1985) 2473. 

[34] A. D. Kennedy, J. Kuti, S. Meyer and B. J. Pendleton, Phys. Rev. D 38 
(1988) 627. 

[35] B. Joo, I. Horvath and K, F. Liu, Phys. Rev. D 67 (2003) 074505 
[arXiv:hep-lat/0112033]. 



18 



[36] I. Montvay and E. Scholz, Phys. Lett. B 623 (2005) 73 [arXiv:hep- 
lat/0506006]. 

[37] M. Weingart, 'Stochastic Estimator of the s Quark Determinant in Full 
QCD Simulation', Diploma work, Uni. Bern, 2007. 

[38] R. Sommer, Nucl. Phys. B 411 839 (1994) [arXiv:hcp-lat/9310022]. 

[39] A. Hasenfratz, P. Hascnfratz, D. Hicrl, F. Niedermayer and A. Schafer, PoS 
LAT2006 (2006) 178 [arXiv:hep-lat/0610096]. 

[40] P. Hasenfratz,V. Laliena and F. Niedermayer, Phys. Lett. B 427 (1998) 
125 [arXiv:hep-lat/9801021]. 

[41] H. Neuberger, Phys. Rev. Lett, 60 (1988) 880; Nucl. Phys. B 300 [FS22] 
(1988) 180. 

[42] G. Martinclli, C. Pittori, C. T. Sachrajda, M. Testa and A. Vladikas, Nucl. 
Phys. B 445 (1995) 81 [arXiv:hcp-lat/9411010]. 

[43] V. Maillart and F. Niedermayer, arXiv:0807.0030 [hep-lat]. 

[44] E. Franco and V. Lubicz, Nucl. Phys. B 531 (1998) 641 [arXiv:hep- 
ph/9803491]. 

[45] K. G. Chctyrkin and A. Rctcy, Nucl. Phys. B 583 (2000) 3 [arXiv:hep- 
ph/9910332]. 

[46] K. G. Chctyrkin, Phys. Lett. B 404 (1997) 161 [arXiv:hcp-ph/9703278]. 

[47] J. A. M. Vermaseren, S. A. Larin and T. van Ritbergen, Phys. Lett. B 405 
(1997) 327 [arXiv:hcp-ph/9703284]. 

[48] T. van Ritbergen, J. A. M. Vermaseren and S. A. Larin, Phys. Lett. B 400 
(1997) 379 [arXiv:hcp-ph/9701390]. 

[49] Y. Aoki et al, Phys. Rev. D 78 (2008) 054510, arXiv:0712.1061 [hep-lat]. 

[50] W. M. Yao et al. [Particle Data Group], J. Phys. G 33 (2006) 1. 



19 



